Black hole growth in the early Universe is 
self-regulated and largely hidden from view 
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Cfche formation of the first massive objects in the infant Universe re- 
^mains impossible to observe dir ectly and yet it sets the stage for the 
'^Subsequent evolution of galaxies-U^l^u While some black holes with 
"glasses >10 9 M Q have been detected in luminous quasars less than one 
pillion years after the Big Bang 3 6S , these individual extreme objects 
' have limited utility in constraining the channels of formation of the 
e arliest black holes. The initial conditions of black hole seed properties 
K^re quickly erased during the growth process 6 . From deep, optimally- 
Q^tacked, archival X-ray observations, we measure the amount of black 
j--4iole growth in 2=6-8 galaxies (0.7-1 billion years after the Big Bang). 
f~^)ur results imply that black holes grow in tandem with their hosts 
(^throughout cosmic history, starting from the earliest times. We find 
i fQiat most copiously accreting black holes at these epochs are buried 
^— ui significant amounts of gas and dust that absorb most radiation ex- 
c ept for the highest energy X-rays. This suggests that black holes grow 
^significantly more than previously thought during these early bursts, 
• and due to obscuration they do not contribute to the re-ionization of 
# pfoe Universe with their ultraviolet emission. 

The Chandra X-ray observatory is sensitive to photons in the energy 
S-aange 0.5-8 keV, which in deep extragalactic observations probes predom- 
inantly accretion onto supermassive black holes 7 . Rapidly growing black 
holes are known to be surrounded by an obscuring medium, which can 
block most of the optical, ultra-violet and even soft X-ray photons^. With 
increasing redshift, at the earliest epochs, the photons observed by Chan- 
dra are emitted at intrinsically higher energies, and therefore less affected 
by such absorption. Current X-ray observations have not been able to indi- 
vidually detect most of the first black hole growth events at z>6 (first 950 
million years after the Big Bang) thus far, except for the most luminous 
quasars 9 at Lx>3x 10 44 erg s~ 4 . While deep X-ray surveys do not cover 
enough volume at high redshift, current wide-area studies are simply not 
deep enough. Hence, the only way to obtain a detectable signal from more 
typical growing black holes is by adding the X-ray emission from a large 
number of sources at these redshifts, which we pursue here. 

We start by studying the collective X-ray emission from the most dis- 
tant galaxies known, at detected by the Wide 
Field Camera aboard the Hubble Space Telescope. These galaxies are as 
massive as today's galaxies ( 10 9 — lO n M0 stellar mass^), and they 



are thus likely to harbor substantial nuclear black holes. None of the 
z>6 galaxies studied in this work are individually detected in the Chan- 
dra X-ray observations. However, we detect significant signals from a 
stack of 197 galaxies at z~6 in both the soft (0.5-2.0 keV; correspond- 
ing to 3.5-14 keV in the rest frame) and hard (2-8 keV; rest-frame 14- 
56 keV) X-ray bands independently. The detection in the soft band is 
significant at the 5-cr level and implies an average observed-frame lumi- 
nosity of 9.2xl0 41 erg s _1 , while in the hard band the stacked 6.8-ct 
signal corresponds to an average luminosity of 8.4 xlO 42 erg s _1 . For 
the sample of galaxies at z~7 we obtain 3-a upper limits for the av- 
erage luminosity in the observed-frame soft and hard X-ray bands of 
4xl0 42 erg s _1 and 2.9xl0 43 erg s 1 respectively. Combining the z~7 
and 2~8 samples the corresponding 3-a upper limits are 3.1 x 10 42 erg s _1 
and 2.2x 10 43 erg s~ 4 in the observed-frame soft and hard X-ray bands. 

A large difference, of a factor of ~9, is found between the stacked 
fluxes in the soft and hard X-ray bands at z ~ 6. This requires large 
amounts of obscuring material with high columns (Nh>1-6x 10 24 cm -2 ) 
to be present in a very high fraction of the accreting black holes in these 
galaxies, in order to explain the large deficit of soft X-ray photons. Since 
this signal derives from the entire population, these results require that al- 
most all sources are significantly obscured. This in turn implies that these 
growing black holes are obscured along most lines of sight, as observed 
in a small subset of nearby objects 39 as well. Such high fraction of ob- 
scured sources at low luminosities is also observed at low redshifts^. This 
large amount of obscuration along all directions absorbs virtually all ultra- 
violet photons from growing black holes. Thus, regardless of the amount 
of accretion in these sources, these active galaxies cannot contribute to the 
early re-ionization of the Universe. Alternatively, it cannot be claimed that 
rapid and efficient supermassive black hole growth in the high-z Universe 
is implausible on the basis of any re-ionization constraints-^. If most of 
the high-redshift black hole growth is indeed obscured as suggested by our 
work, several current constraints on the lifetime and duty cycle of high-z 
accreting black holes need to be revisited and revised. 

Assuming that the X-ray emission is due to accretion onto the central 
black hole, the accreted black hole mass density can be directly derived 
from the observed X-ray luminosity, as described in the supplementary 
information. Extrapolations of Active Galactic Nuclei (AGN) luminosity 
functions 58 measured at significantly lower redshifts, z<3, are consistent 
with the observed accreted black hole mass density at z>6, as can be seen 
in Figure 1. This directly leads to two further conclusions: the space den- 
sity of low luminosity sources, Lx<l0 44 erg s _1 , does not evolve signif- 
icantly from z~l to z~6-8, i.e. over more than 5 billion years. Second, 
at higher luminosities, the extrapolation of lower-redshift AGN luminosity 
functions leads to an overestimate of the observed source density in opti- 
cal surveys^. This discrepancy can be resolved if the shape of the AGN 
luminosity function evolves strongly in the sense that there are relatively 
fewer high-luminosity AGN at z>6 in comparison to the z<3 population. 
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Another possibility is that the number of obscured sources, relative to un- 
obscured quasars, increases with redshift, such that most of the highly ob- 
scured systems are systematically missed in these optical studies. This is 
strongly supported by observations of quasars at lower redshifts, z<$^. 
We cannot rule out either of these scenarios at present due to the rel- 
atively small cosmological volume studied, in which the extremely rare 
high-luminosity AGN are absent. 

Our measurements and upper limits for the accreted black hole mass 
density up to z~8.5 (~600 million years after the Big Bang) constrain the 
nature of black hole growth in the early Universe. Two critical issues for 
AGN and the supermassive black holes powering them are how the first 
black holes formed, and how they subsequently grew accreting mass while 
shining as AGN. The strong local correlation between black hole mass and 
galaxy bulge mass observed at 2~C 20 21 , is interpreted as evidence for self- 
regulated black hole growth and galaxy-black hole co-evolution ' 22 . This is 
currently the default assumption for most galaxy formation and evolution 
modelW 

The origin of the initial "seed" black holes remains an unsolved prob- 
lem at present. Two channels to form these seeds have been proposed: 
compact remnants of the first stars, the so-called population III stars 2 ^, 
which generate seeds with masses ~ 10-1,000 M© and from the direct 
gravitational collapse of gas-rich pre-galactic disks, which leads to sig- 
nificantly more massive seeds with masses in the range M~10 5 MpPlZl. 
By construction, the masses of seeds that form from direct collapse are 
correlated to properties of the dark matter halo and hence properties of the 
galaxy that will assemble subsequently. 

To interpret our finding, we explore a theoretical framework for the 
cosmic evolution of supermassive black holes in a ACDM cosmology. We 
follow the formation and evolution of black holes through dedicated Monte 
Carlo merger tree simulations. Each model is constructed by tracing the 
merger hierarchy of dark matter halos in the mass range 10 11 — 1O 15 M0 
backwards to z = 20, using an extended Press & Schechter algorithm^! 
The halos are then seeded with black holes and their evolution is tracked 
forward to the present time. Following a major merger (defined as a merger 
between two halos with mass ratio > 0.1), supermassive black holes ac- 
crete efficiently an amount of mass that is set by a "self-regulated" model 
(where the accreted mass scales with the fourth power of the host halo 
circular velocity and is normalized to reproduce the observed local corre- 
lation between supermassive black hole mass and velocity dispersion) or a 
"un-regulated" model, where the supermassive black hole simply doubles 
in mass at each accretion episode. See the Supplementary Information for 
additional details. 

Our observational results provide strong support for the existence of a 
correlation between supermassive black holes and their hosts out to the 
highest redshifts. In Figure 1, we compare both unregulated and self- 
regulated black hole growth models with our observations, and find that 
physically motivated self-regulation growth models are highly favored at 
all redshifts, even in the very early Universe. Un-regulated models (for 
instance wherein black holes just double in mass at each major merger) 
are strongly disfavored by the data. This indicates that even in the first 
episodes of black hole growth there is a fundamental link between galaxy 
and black hole mass assembly. 

As shown in Figure 1, once a standard prescription for self-regulation 
(as described before) is incorporated, both seed models are consistent with 
our current high-z observations. Detection of an unbiased population of 
sources at these early epochs is the one metric that we have in the fore- 
seeable future to distinguish between these two scenarios for the origin 
of supermassive black holes in the Universe. In Figure 2, we present the 
predicted cumulative source counts at z>6 for the models studied here. 
Based on these models, ultra-deep X-ray and near-infrared surveys cov- 
ering at least ~ 1 deg 2 are required to constrain the formation of the first 
black hole seeds. This will likely require the use of the next generation of 
space-based observatories such as the James Webb Space Telescope and 



the International X-ray Observatory. 
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Figure 1 | Accreted black hole mass density as a function of redshift. The 
gray rectangle shows the range of values allowed by observations of z~0 
galaxieJ^l. The data points at z~2 correspond to the values obtained 
from Chandra observations of X-ray detected AGN and luminous infrared 
galaxies 54 , while the measurement at z~6 and the upper limits at z=7-9 
show the results described in this work (red and black data points from the 
observed-frame soft and hard X-ray band observations respectively). Vertical 
error bars represent 1 s.d. while the horizontal ones show the bin size. The 
black solid line shows the evolution of the accreted black hole mass density 
inferred from the extrapolation of AGN luminosity functions measured at lower 
redshifts 58 . We over-plot the predictions of black hole and galaxy evolution 
models 30 for non-regulated growth of Population-Ill star remnants (cyan line) 
and direct-collapse seeds (green). The red and blue lines show the predicted 
BH mass density if self-regulation is incorporated. 
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Figure 2 | Cumulative number of sources as a function of redshift for individual 
X-ray detections. This calculation assumes the X-ray flux limit of the 4 Msec 
CDF-S Chandra observations. The horizontal dotted line shows the number 
density required to individually detect one source in the area considered in 
this work at z >7. Models are described in the supplemental material and 
labelled in the figure (Pop III, / E dd = 1; D.C., / E dd = 0.3; D.C., / E dd = 0.3, 
x2). Note: model Pop III, / E dd = 1, x2 has no detectable source. To 
distinguish between these models for early black hole formation will require a 
deep multiwavelength survey covering at least ~1 deg 2 . 
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The Observations 

The main observational data samples for this work are the 371 galaxy 
candidates at z~6 selected using the optical and near-IR Lyman break 
technique 31 , together with 66 z^l and 47 2-~8 galaxy candidates^! a rj 
of them in the Chandra Deep Field South (CDF-S) field. We then com- 
plemented this sample with the 151 z~6 candidates in the Chandra Deep 
Field North (CDF-NjM The accuracy of this drop-out technique selec- 
tion has been recently confirmed by spectroscopic observations of some of 
these sourcePEH. The contamination by foreground sources in these sam- 
ples appears to be very small, ~109wSE3. We then used the 4 Msec Chan- 
dra observations of the CDF- S^and the 2 Msec Chandra data available on 
the CDF-N^in order to search for signatures of supermassive black hole 
(SMBH) accretion in these sources. None of these galaxies are detected 
individually in the X-ray data. However, Chandra data are uniquely suited 
to perform stacking, which allows us to reach much fainter flux limits. In 
order to maximize the signal-to-noise of these measurements we used an 
optimized X-ray stacking scheme, described in detail in Appendix A. We 
stack independently in the soft, 0.5-2 keV, and hard, 2-8 keV, observed- 
frame Chandra bands. At z~7, they correspond to rest-frame energies of 
4-16 and 16-64 keV respectively. While at these high rest-frame energies 
the effects of Compton-thin obscuration (Nh <10 24 cm~ 2 ) are negligible, 
we have to consider the possibility of higher absorption column densi- 
ties. We restricted our stack to sources closer than ~9' from the average 
aim point of the Chandra observations in order to have an optimal extrac- 
tion radius smaller than 7" (Fig. S|3|. We further removed sources with 
an X-ray detection closer than 22" to avoid possible contamination in the 
background determination. 

In the sample of candidates at z~6 we stacked a total of 197 sources, 
151 in the CDF-S and 46 in the CDF-N. This corresponds to a total 
exposure time of ~7xl0 8 seconds (~23 years). We found significant 
detections, >5-er, in both the soft and hard bands independently. In 
the soft band we computed a count rate of 3.4±0.68x 10~ 7 counts s _1 
per source, which corresponds to an observed-frame soft band flux of 
2.3x 10~ 18 erg cm _2 s _1 . Converting it to the rest-frame hard band we 
obtain a flux of 1.9xl0 -18 erg cm _2 s~ 1 . In the observed-frame hard 
band we measure a count rate of 8.8±1.3x 10 -7 counts s~\ which corre- 
sponds to a 6.8-er detection. The average flux in the observed-frame hard 
band is2.1xl0~ 17 erg cm _2 s _1 , which converted to the rest-frame hard- 
band corresponds to 1.7x 10 -17 erg cm _2 s -1 , Stacked Chandra images 
for this sample in the soft and hard X-ray bands are shown in Fig. S|4] 

There is a factor ~9 difference between the fluxes measured in the 
observed-frame soft and hard bands. Assuming a power-law X-ray spec- 
trum with an intrinsic spectral slope F=1.9, typical of AGN 38 , the expected 
flux ratio between the observed-frame soft and hard bands for an unob- 
scured source is expected to be ~1.7 (Fig. Sj5}. The only explanation for 
the relatively large flux ratio in the hard to soft bands is very high levels 
of obscuration. As can be seen in Fig. S|5] at 2~6 a minimum column 
density of N_h~10 24 cm -2 , i.e. Compton-thick obscuration, is required. 
Given that this ratio is observed in the stack, this implies that there are 
very few sources with significantly lower levels of obscuration, which in 
turn means that these sources must be nearly Compton-thick along most 
directions (~4-7r obscuration). Similar sources have also been observed 
in the local Universe 39 but are likely rare. Furthermore, up to z~3 it has 
been shown before 40 42 that the fraction of obscu red AGN increases with 
decreasing luminosity and increasing redshifl 43 44 . Hence, it is not entirely 
surprising that the sources studied here, given their low luminosities and 
high redshifts, are heavily obscured. In fact, the discovery of a Compton- 
thick AGN at z~5 selected using the drop-out technique has been recently 
reported 4 ^! 

The corresponding average rest-frame 2-10 keV luminosity, derived 
from the observed-frame hard band, is 6.8xl0 42 erg s _1 . Since none 



10 



-i i i | i i i | i r- 

Soft band (0.5-2 keV) 
Hard band (2-8 keV) 



■o 

CO 



ft 
o 




2 4 6 

Off-axis angle [arcmin] 



10 



Figure 3 | Optimal extraction radius as a function of off-axis angle, measured 
from the average Chandra pointing center. As described in the appendix A, 
these lines were measured by optimizing the function f(8, r)/r. A minimum 
radius of 1 " was assumed in order to avoid flux losses due to astrometric 
problems and pixel aliasing. Red and blue lines show the optimal radii for the 
soft and hard band respectively. Sources with optimal radii greater than 7" 
(dashed horizontal line) are not considered in the stack, given their low expected 
contribution to the integrated signal-to-noise. 




Figure 4 | Stacked Chandra images for the z=6 galaxy sample in the soft (left 
panel) and hard (right panel) X-ray bands. The detections are significant at 
the 5 and 6.8 -a levels respectively. Each image is 30"x30". The white 
circle at the center of each image has a radius of 3". Images were adaptively 
smoothed using a minimum scale of 3 pixels, a maximum scale of 5 pixels and 
minimum and maximum significances of 3 and 6 respectively. 
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Figure 5 | Expected ratio of the observed-frame hard to soft flux as a function 
of obscuring neutral Hydrogen column density (N H ). The black solid line was 
derived assuming an intrinsic power-law spectrum with slope r=1.9 and pho- 
toelectric absorption. The gray zone shows the measured ratio for the stack of 
galaxies at z~6 and the ±1 s.d. limits. A column density of N^c^lO 24 cm~ 2 , 
i.e. Compton-thick obscuration, is required to explain the observed hard to soft 
X-ray flux ratio. 

of these sources are individually detected in X-rays, we conclude that 
at least 30% of the galaxies in this sample contain an AGN. Mul- 
tiplying by the number of sources we obtain a total luminosity of 
1.34xl0 45 erg s _1 . The total area surveyed is ~310 arcmin 2 , or 0.086 
deg 2 . Hence, the integrated AGN emissivity at z~6 derived from this 
sample is 1.6x 10 46 erg s _1 deg -2 . 

Similarly, we stacked a total of 57 sources by combining the sam- 
ples at z>l and z>8. For the galaxies in this redshift range, we 
found no significant detection in the observed-frame soft or hard bands. 
The 3-cr upper limits are 6.9xl0 -7 cts s _1 and 1.4xl0 -6 cts s _1 , 
which corresponds to flux upper limits of 4.6x 10~ 18 erg cm" 2 s" 1 and 
3.3xl0~ 17 erg cm~ 2 s _1 respectively. Assuming an average redshift 
2=7.5 and converting to the rest-frame hard band, 2-10 keV, we obtain 
average luminosities of 3.1 x 10 42 erg s _1 and 2.2x 10 43 erg s -1 from the 
observed-frame soft and hard band respectively. Because we only have up- 
per limits in both bands, we cannot constrain the presence of obscuration 
in this sample, however assuming that these AGN are as heavily obscured 
as the z~6 sample, the shallower limit from the hard band is actually the 
more stringent constraint. The area covered in the z>l observations is 40 
arcmin 2 , or 0.011 deg 2 . Hence, the upper limits to the integrated AGN 
emissivity are 5.2xl0 45 erg s _1 deg -2 and 4.3xl0 46 erg s _1 deg -2 re- 
spectively. 

Finally, we also stacked the galaxies in the z~7 sample separately. In 
this case we obtained 3-cr upper limits of 6.9xl0 -18 erg cm _2 s _1 and 
5. Ox 10 -17 erg cm _2 s _1 in the observed frame soft and hard bands re- 
spectively. These correspond to average observed-frame luminosities of 
4. Ox 10 42 erg s" 1 and 2.9x 10 43 erg s _1 . Dividing by the survey area and 
multiplying by the number of sources we obtain 3-cr upper limits for the in- 
tegrated AGN emissivities from these galaxies of 4.3 x 10 45 erg s -1 deg -2 
and 2.9 x 10 46 erg s" 1 respectively. 



Black Hole Mass Density Determination 

We compute the observed integrated black hole mass density using an 
updated v ersion of the "Soltan'' 4 ^ argument. Following the standard 
derivationJ2E3 W e have that the integrated black hole mass density is given 
by: 

r oo it /*°° 1 £ 

Pbh(z)= — dz / -^ r L bo i9(L,z)dL, (1) 
J z dz J ec 

where 

Lboi — k corr Lx (2) 

and k corr is the bolometric correction for the rest-frame hard X-ray 
band. In order to compute pbh(z) we need to make some assumptions, 
since the AGN luminosity function (LF) at z>6 is unknown, and only the 
integrated luminosity is known. First, we assume that at z>6 the AGN 
LF does not depend strongly on redshift, i.e., ^(L, z)~^(L). We further 
assume that the efficiency and bolometric corrections are constant (i.e., 
independent of luminosity). While we know that the bolometric correction 
depends on luminosity—^H our sample most likely does not span a wide 
luminosity range, and hence this assumption is reasonable. Therefore, we 
have, 

Pbh(z) = (1 ~ £) 2 fcc ° rr r ^dz r L x *(L)dL. (3) 
ec J z dz J 

The second integral on the right hand side can be determined from 
the observed integrated AGN emissivity. We only need to convert from 
a density per unit of sky area (the observed value) to a co-moving AGN 
emissivity per unit of volume. This can be done easily by dividing the 
observed value by the co-moving volume at the redshift range covered by 
each sample. In principle, a correction should be made to account for 
the contribution of high-luminosity sources not present in the relatively 
narrow fields considered in this work. However, as we will show below, 
this correction should be small, ~l-2%, and thus can be safely ignored. 

The bolometric correction fc corr has been estimated by several au- 
thors in the past. For high-luminosity sources (quasars), a value of 35 
was measureo 50 . For low-luminosity sources, L^IO 43 erg s _1 , values 
of ~ 10-20, were estimated^, while newer calculation^ report a value of 
~25 for sources with L x ~10 42 erg s _1 and ~40 at L x ~10 43 erg s" 1 . 
Given the low luminosity of the sources in our sample we assume a value 
of fc corr =25, with an estimated uncertainty of a factor of ~2. The main 
factor contributing to this uncertainty is the determination of the contri- 
bution of the X-ray emission relative to the ultraviolet (where most of the 
bolometric output is found). Observational studies show that this factor 
depends on luminosity but remains constant with redshift up to 2~(P2I 

Assuming a constant radiation efficiency e=0.1 we obtain from the 
combined 2=7-8 sample 3-cr upper limits of pbh<T0S MqMpc -3 
(possibly affected by obscuration) and p.BH<5883M0Mpc~ 3 from the 
observed-frame soft and hard band respectively. 

For the sample of galaxies at z~7 only, assuming that all these galaxies 
lie between zi=7 and 22=8, we obtain that psir<1117 Mq Mpc -3 from 
the observed-frame soft-band and pBH <7595 Mq Mpc -3 from the hard 
band. Finally, for our sample of galaxy candidates we obtain from the 
stacked detection in the observed-frame hard band that psH=5005±751 
Mq Mpc" 3 (l-cr). 

In Figure S|6] we plot the integrated accreted black hole mass density 
as a function of redshift. The upper limits at z~7-8 and measurement at 
z^6 obtained from our work are shown together with the observations in 
the local Universe^ and the values derived from Chandra observations of 
X-ray detected AGN and luminous infrared galaxies at 

Accretion Models: key features 

We investigate the formation and evolution of black holes via cosmolog- 
ical realizations of the merger hierarchy of dark matter halos from early 
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Figure 6 | Accreted black hole mass density as a function of redshift. The 
gray rectangle shows the range of values allowed by observations of z~0 
galaxies^! The data points at 2~2 correspond to the values obtained 
from Chandra observations of X-ray detected AGN and luminous infrared 
galaxies 54 , while the measurements at z~6 and the upper limits at z=7-9 
show the results described in this work (red and black data points from the 
observed-frame soft and hard X-ray band observations respectively). T he red 
and blue lines show the values inferred from AGN luminosity functions^ 5 -! 5 -^ 
while the point at z=6-7 was obtained from the quasar luminosity function 57 . 
The black line assumes the hard X-ray AGN LF of Ueda et alP, as modified 
by Treister et al! 5 -^. 

times to the present in a ACDM cosmology. The main features of the 
by now fairly standard SMBH evolution models that a re used to interpret 
the data have been discussed in detail elsewher e 1 59 * 6 ' !. We briefly outline 
some of the key features here. In this work, two "seed" formation mod- 
els are considered : those deriving from population-Ill star remnants (Pop 
III), and from direct collapse models (D.C.). The main difference between 
these two models lies in the mass function of seeds. Pop III seeds are 
light weight (few hundred solar masses) and form abundantly and early 
(roughly one per comoving cubic Mpc at z ~ 20). D.C. seeds are more 
massive (10 4 — 10 6 solar masses), but rarer (a peak density of 0.1 per co- 
moving cubic Mpc at z ~ 12). We summarize below the relevant and 
standard assumptions that go into our modeling of SMBH growth. 

Essentially in this scheme central SMBHs hosted in galaxies accumu- 
late mass via accretion episodes that are triggered by galaxy mergers. Ac- 
cretion proceeds in one of two modes: self-regulated or un-regulated. For 
each SMBH in our models we know its mass at the time when the merger 
starts (Afi n ), and we set the final mass through the self-regulated or un- 
regulated prescription. These two models differ by the amount of mass a 
SMBH accretes during a given accretion phase. In the context of the cur- 
rently supported paradigm for structure formation, growth of structure in 
the Universe occurs hierarchically and via copious merging activity. Our 
models rely on the following assumptions: 

• SMBHs in galaxies undergoing a major merger (i.e., having a mass 
ratio >1:10) accrete mass and become active. 

• In the self-regulated model, each SMBH accretes an amount of 
mass, corresponding to 90% of the mass predicted by the local A/bh — & 



relation 62 , 

M fln = M iB + 0.9 x 1.3 x 10* {^^y 24 M ; (4) 

the 90% normalization was chosen to take into account the contribution 
of mergers, without largely exceeding the mass given by the A/bh — a 
relation for SMBHs at z — 0. Here a is the velocity dispersion of the host 
after the merger. We adopt a — V c /v3, where V c is the virial velocity of 
the host dark matter halo 63 ' 64 ! . A SMBH is assumed to stop accreting once 
it reaches the value given by the Mbh — cr relation. 

•In the unregulated mode (x2) we simply set A/fi n = 2 x M- m , that 
is, we double the mass during each accretion episode. 

• The rate at which mass is accreted scales with the Eddington rate 
(/Edd) for the SMBH, where /Edd is the accretion rate in units of the 
Eddington rate. We adopt /Edd = 0.3 for D.C. models and /Edd = 1 for 
Pop III models in order to reproduce the mass density at z — 0. As Pop III 
seeds are lighter, they need to accrete more mass over the course of their 
cosmic history to become supermassive. 

• The lifetime of an AGN depends on how much mass it accretes dur- 
ing each episode. Given the initial mass of a SMBH, M; n , and the amount 
of mass it accretes, we can calculate its final mass, Mfi n , at the end of the 
active phase (Equation 4). For a given Eddington fraction, /Edd, the mass 
of the SMBH grows with time as: 

A/fin = M- m exp ( / Edd — € ) (5) 

\ £ lEdd / 

where e is the radiative efficiency (e ~ 0.1), tEdd = A/BHC 2 /I/Edd = 
47r°Gm ~ 0.45 Gyr (c is the speed of light, or is the Thomson cross 
section, m p is the proton mass). Given Mi n and Mjj n , then a SMBH 
is active for a time £agn = jt^ ln(M fln /M in ). Note that the un- 
regulated prescription (a SMBH mass doubles at each accretion episode) 
corresponds to assuming that all SMBHs shine as AGN for the same time 
at a fixed accretion rate (e.g., 100 Myr for an accretion rate /sdd = 0.3). 

In summary, we study and compare two self-regulated models (Pop 
III, /Edd = 1; D.C., /Edd = 0.3) and two unregulated models (Pop III, 
/Edd = 1, x2; D.C, /Edd = 0.3, x2). We use these models to derive 
the mass density in black holes in a cosmic volume, by summing over all 
existing black holes at a given redshift and normalizing by the comoving 
volume. We also calculate number counts (Figure 2 in the main article) 
imposing selection criteria that match the observations we have analyzed. 
The limiting luminosity is 4.2 x 10 41 erg s _1 for the stacked sample, i.e., 
if all the galaxies are emitting at this limit then we should detect them 
in the stack. For a single source, the more appropriate limit is this value 
multiplied by the number of sources in the stack, which is equivalent to a 
single source contributing all the signal. In that case, the luminosity limit 
is 3.1 x 10 43 erg s" 1 and the flux limit is 4.95 x 10~ 17 erg cm -2 s" 1 . 
In that case, the source will be individually detected. We therefore select 
all accreting black holes in the theoretical sample at z > 6 that meet this 
flux criterion. From the Eddington ratio we calculate the bolometric lumi- 
nosity and then apply a K-correction of 25 to derive the X-ray flux at the 
appropriate redshift. 

Comparison with AGN Luminosity Functions 

We first compare the observational results presented with the expectations 
from integrating existing measurements or extrapolations of the observed 
AGN luminosity functions from earlier work. Integrating the hard X- 
ray LP^ from z=l to 10 we obtain an integrated luminosity density of 
2.3 x 10 46 erg s _1 deg~ 2 , in good agreement with the values obtained from 
the stacking in the hard X-ray band, which are less affected by obscura- 
tion. In contrast, from the LF of z~6 quasars 57 , we expect a steep decline 
in the AGN number density. 

Comparison of observations of the accreted SMBH mass density with 
predictions obtained from extrapolations of existing AGN luminosity func- 
tions to high redshifts and/or low luminosities provides contrasting results. 
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While there is in general good agreement up to z~3, differences of up to 
~2 orders of magnitude are found at z>6. This can be explained as the 
luminosity functions of Hopkins et alP^and Willott et al. incorporate a 
density evolution of the form 10~ 0A7z , which is required to account for 
the observed evolution of the highest luminosity quasars at z~5-6 in the 
SDSS^! A similarly steep evolution from jz~3 to z~5 has been reported, 
based on a compilation of Chandra X-ray observations 5 * . However, it is 
worth pointing out that the accreted SMBH mass density inferred from 
this LF, presented in Figure S[6] is significantly lower than observations 
and other LFs at all redshifts, which may be explained by the relatively 
high spectroscopic incompleteness (~50%) of this sample. Incomplete- 
ness is a particularly important issue and limitation at high redshift, where 
sources are faint in the optical bands. In contrast, the AGN LF of Treister 
et al.^ is in good agreement with the observations at all redshifts, and it 
is therefore the only one plotted in Figure 1 of the main paper. This work 
is based on the Ueda et al. 40 hard X-ray AGN LF modified to incorporate 
an increasing number of obscured sources with increasing redshifl 4 ^ and 
reducing the relative number of Compton-thick sources by a factor of ~4, 
consistent with the observed space density of these sources in the all-sky 
INTEGRAL and Swift surveys. 

Explaining the accreted SMBH density at z>6 inferred from X-ray 
stacking observations would require that the comoving space density of 
low-luminosity AGN, Lx<W 44 erg s - , remains nearly constant up to 
z~3@2| At the same time, the discrepancy with the observed density of 
high-luminosity sources in optical surveys can be resolved if there is a 
strong evolution in the number of high-luminosity sources. For example, a 
decline in the number of quasars with redshift given by l0-° 47z has been 
measured^. In this case, the contribution from high-luminosity sources 
to the integrated accreted black hole mass density will be very small, ~1- 
2%. Alternatively, the observed results could be explained if the relative 
number of heavily-obscured quasars increases strongly with redshift, as 
measured up to z~3^. 

Implications for re-ionization 

Two sources of UV radiation are expec ted to contribute to re-ionization: 
star forming galaxies and quasar J 67 * 69 ! The conventional view is that 
galaxies dominate hydrogen re-ionization occurring at high redshift 70 , 
z ~ 9, and quasars dominate helium re-ionization occurring at a later 
time, z ~ Jp^. This result is based on optically-selected quasar luminosity 
functions 5 ^" that show, as discussed in the previous section, a steep drop 
of the quasar population at z>4^. However, in ab initio models of black 
hole evolution through cosmic history^ the amount of black hole growth 
required to explain the bright end of the luminosity function of quasars 
at z — 3 — 6 implies a substantial contribution of quasars to hydrogen 
re-ionization, provided that UV and X-ray radiation can escape into the 
intergalactic medium. Our findings resolve this tension, as regardless of 
the amount of accretion occurring onto black holes, the bulk of the emis- 
sion from the population is obscured, and UV photons cannot escape. In 
4-7T Compton thick sources only the highest energy photons can escape. 
We provide a simple estimate here of the contribution of Compton thick 
quasars to re-ionization. From our models we can extract information on 
dp, the accreted mass density on black holes as a function of cosmic time 
(Figure 1 in the main article shows the integral quantity). Let us assume, 
optimistically, that an escape fraction f esc — 0.04 of the bolometric lumi- 
nosity (here we just assume / esc = 1 jk corr ) goes into hard X-ray photons 
with mean energy E 1 = 10 keV that can escape from the quasar host. 
All these photons contribute to primary and secondary ionizations. The 
number of ionizing photons per hydrogen atom can be evaluated from the 
following differential equation: 

dx= J^ d + /si(s) _ 

p H E-y F p H UAeV F t rec ' 



where pu is the hydrogen cosmic density^, and the recombination 
timescale for hydrogen is t rec ~ 0.3[(1 + z)/4] -3 GyM Here the first 
term of the right hand side of the equation includes primary ionizations by 
ionizing photons, the seco nd term accounts for secondary ionizations from 
energetic ionizing photons I22E4I ant ] t ne third term accounts for recombi- 
nations. Here f S i(x) w 0.35(1 - a; ' 4 ) 1 ' 8 - 1.77 )°' 4 x 02 (l - 

x 0A ) 2 , as long as x < 1, and x is less than unity in all our models down 
to z ~ 1. 

The number of ionizing photons per hydrogen atom is shown in the 
accompanying Figure ^7] where we see that models that successfully re- 
produce our observations fail to reionize hydrogen by 2-3 orders of mag- 
nitude at z >6, notwithstanding significant black hole growth. Without 
taking obscuration into account, these accreting black holes would con- 
tribute substantially to re-ionization, producing one ionizing photon per 
atom by z = 8^. 

However, we stress here that if quasars do not contribute to hydro- 
gen re-ionization, it is not because quasar activity drops steeply at z>4 
as previously suggesteG 71 , but instead due to the presence of obscuring 
material that absorbs UV and soft X-ray radiation. As seen in Figure ^7] 
the models that are consistent with the observations at z>6 and include 
self-regulation under-produce ionizing photons at these epochs. Hence, re- 
ionization of the early Universe is most likely due to star forming galaxies 
and not growing black holes. However, this is still vigourously debated. 
As previously showrP 5 -', observed z~7 galaxies cannot provide enough 
hydrogen-ionizing photons unless some of the galaxy properties, such as 
escape fraction, metallicity, initial stellar mass function or dust extinction, 
evolve significantly from z~7 to the local Universe. Similarly, the extrap- 
olation of the galaxy luminosity function at the faint end plays a major role 
in whether these sources can re-ionize the Universe^. 
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Figure 7 | Number of ionizing photons per hydrogen atom. Models are labelled 
in the figure (Pop III, / E dd = 1; Pop III, / E dd = 1, x2; D.C., / E dd = 0.3; 
D.C., / E dd = 0.3, x2). A ratio of at least one ionizing photon per hydrogen 
atom is required to re-ionize the Universe. As can be seen in this figure, most 
models predict a ratio for the early growing black holes that is lower than this 
value by 2-3 orders of magnitude at z>6. 
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Appendix A: Optimized X-ray stacking 

We start by assuming that all sources are equally bright in the X-rays, with 
count rate R 3 , Rh in cts/Ms in the soft and hard bands, respectively. The 
background levels in Chandra for the ACIS-I detectors as measured in 
the CDF-S 4 Msec data are B s =0.049 cts/Ms/pixel in the soft band and 
Bh =0.160 cts/Ms/pixel in the hard one, in good agreement with the val- 
ues of 0.066 cts/Ms/pixel (soft), and 0.167 cts/Ms/pixel (hard) measured 
in the CDF-S 2 Msec data^l The background level was found to vary by 
~5% (standard deviation) across the field. To account for this, the back- 
ground value was determined for each source independently by measuring 
the average count rate on a annulus around the source position with inner 
radius 2r, where r is the aperture radius in pixels and outer radius fixed 
to 22". A 3-cr clipping was used for the background determination. Then, 
an aperture of radius r pixels has area ivr 2 and hence background counts 
B k ivr 2 < E(r) > k where < E(r) > k stands for the average of the expo- 
sure time Eij over all pixels inside this aperture around the kth source and 
B k is the background level for that source. The signal contained within this 
aperture can be found by interpolating the enclosed energy profiles, which 
give us the values of r corresponding to various fractions / of R as a func- 
tion of offaxis angle 8, e.g. r(8, f — 0.95). We interpolate between these 
curves to determine the desired function, f(8,r), producing a expected 
signal for the fcth source of S k = Rf(8k, r) < E(r) > k . We assume that 
the noise receives Poisson contributions from both the background and 
the object counts, N k = ^/(Rf(O k , r) + Birr 2 ) < E s (r) >k- Hence the 
signal-to-noise is given by 

_ Rf(8 k ,r)y/<E(r) > k 



N J k ^/Rf(8 k ,r) + Bnri 

For each source, we choose the value of r that maximizes this func- 
tion. In our case, each individual source is dim enough that background 
fluctuations dominate the noise, and assuming that Eij is slowly varying, 
the function to be maximized is just f(8 k ,r) /r. The resulting values of 
rasa function of off-axis angle for the soft and hard bands are shown in 
Figure Sj3] A minimum radius of 1" was assumed in order to avoid flux 
losses due to astrometric problems and pixel aliasing. 

We will stack the estimated source count rates, R k , rather than the 
observed source counts, C k , since the latter depends on the PSF and expo- 
sure time at each source position but the former does not. This is superior 
to performing photometry on the stacked image, which requires using a 
constant aperture for all sources. The standard technique in the literature 
of performing photometry on the stacked image and dividing by the num- 
ber of sources is both sub-optimal and biased; the bias could be fixed by 
instead dividing by the effective number of sources given the PSF and ex- 
posure information. 

While a straightforward averaging ("stack") of the source count rates is 
unbiased, it is still not optimal because it gives equal weight to each count 
rate even though some were measured with lower SNR. We can derive 
the weights that optimize S/N of the combined stack by noting that the 
combined stack will have S = Yli Wi ^i an d N 2 = J2 i w 2 N 2 (since 
the weighted sum is just an unweighted sum of new objects with signal 
WiSi and noise WiNp^). It is easy to show that the optimal weight is 
then Wi = Si /N 2 which has the required behavior that multiplying all the 
weights by a constant k preserves the final S/N. 

Then, we obtain that: 

s = Ej/wjSj = EjgM = /v s 2 /n 2 m 

N vwn? vE.sf/N 2 d 1 

Hence the optimal weights cause S/N to add in quadrature so it will never 
formally decrease. However, in practice there is little benefit to including 
objects that offer individual S/N of less than 10% of the typical objects 
in the field, which is why we only stacked sources closer than 9' from the 
average Chandra aim point. 
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